Questions

This project set out to measure the direct and indirect effects of spatiotemporal variation in salinity on intertidal community structure in the Strait of Georgia, British Columbia. We asked the following questions:

How does the intertidal community vary in low vs. high salinity sites across the span of a summer?
What is the relative impact of grazers in structuring rocky intertidal systems in low vs. high salinity sites?
What is the salinity tolerance of limpets? Is this impacted by emersion from hyposaline conditions during low tides? or by local adaption to salinity regimes?
What is the salinity tolerance of Ulva?

Hypotheses

We predicted that observed differences in community structure between consistently high salinity environments and periodically low salinity environments would be driven by trophic effects as well as physiological stress. We hypothesized that limpets and other marine gastropods are disproportionately affected by the annual decrease in surface salinity experienced in the coastal areas near the mouth of the Fraser River. We further predicted that the abundance of gastropods in periodically low salinity environments would be lower than that of consistently high salinity areas around the Southern Gulf Islands. This reduction in gastropod grazers would lead to a reduction in grazing pressure on palatable green algae, biofilms and algal spores, resulting in greater algal cover in low salinity environments. Increased algal abundance has been shown to negatively impact barnacle species (Farrell 1991), and we therefore predicted lower barnacle abundance in low salinity sites. To test these predictions, we combined laboratory salinity tolerance trials on limpets and green algae with transect surveys and limpet exclusion/inclusion experiments at three sites within each of the two salinity regions.

Methods

Site description

We conducted field studies at three sites within each of two regions with contrasting salinity regimes: West Vancouver, which experiences reduced salinities during the summer, and the Southern Gulf Islands, which experience consistently high salinities year-round. West Vancouver is located approximately 30km north of the Fraser River outflow, and the Southern Gulf Islands are located approximately 30 km southwest (Fig. 1). Gulf island sites are located on the southwest side of the island chain, and are not exposed to the Fraser River plume. The two areas are similar in terms of climate and topography. Sea surface temperature in the two regions is comparable, ranging from 5.0 to 18.5 in West Vancouver and 6.0 to 18.5 in the Gulf Islands (Fisheries and Oceans Canada 2009). The tidal range is greater in West Vancouver, with extreme high tides reaching 4.7 m above Canadian chart datum (approximated as the lowest astronomical tide), compared to 3.4 m above chart datum in the Gulf Islands. All sites were composed of granite rock except HS1, which was sandstone. The slope of the rock face ranged from 1 to 37º and aspect ranged from 40 to 320º east of magnetic north.

Figure 1: Map of the study region. Low salinity sites are located in West Vancouver and high salinity sites are located in the Gulf Islands.

Experimental design

Transect Surveys
Surveys were conducted once per month during low tide from May to August 2011, at each of the six study sites. Because the tidal range differs between the two areas, surveys were conducted at the vertical height corresponding to approximately 30% submersion time. This occurs at 2.1 m in the Southern Gulf Islands and 3.0 m in West Vancouver. Ten meters of transect tape were laid across the selected area and eight randomly selected points were surveyed using a 25x25 cm quadrat. Mobile invertebrates were counted and sessile invertebrates and algae were measured by percent cover.

Field Exclusion Experiment
Seven sub-sites within each of the six study sites (HS1, HS2, HS3, LS1, LS2 and LS3) were manually cleared of organisms. Within each sub-site, a limpet inclusion, exclusion and control plot were set. Inclusions and exclusions were formed by securing two copper fences, 2.5 cm high and 28.5 cm in diameter, to the rock face using Quickcrete® quick drying cement. Copper enclosures/exclosures of this type have been shown to be effective barriers to limpets (Harley 2002) and partial barriers to periwinkles (Harley 2006). Four Lottia pelta, 20±5 mm in diameter, were collected from the nearby shoreline within each site and placed into rings designated as limpet inclusion treatments. This density of 0.63 limpets per 100 cm², approximately corresponds to the average density of limpets in low salinity sites as determined by previous survey data (0.55 per 100 cm² in low salinity sites and 3.65 per 100 cm² in high salinity sites). L. pelta were determined to be the largest grazer present in the mid-intertidal in both salinity regions, and body size of grazers has been shown to positively correlate with grazing pressure (Geller 1991). Any other grazers found inside rings were removed. Limpet treatment densities were maintained every two weeks by adding or removing limpets as necessary. One circular plot within each area, also 28.5 cm in diameter, was marked with steel bolts and served as a control. The position of control and treatments was randomized within each subsite. Copper controls were not included in this study, as previous work has shown that partial copper treatments lead to partial effects which are difficult to interpret (Johnson 1992).

Sampling occurred once per month during low tide from May to August. A 10x10 cm quadrat was used to count mobile invertebrates and barnacles and estimate percent cover of algae and mussels within each treatment. Salinity samples were taken at each sampling event and measured using a refractometer (S/Mill-E, Atago Inc.).

Statistical Analyses

All analyses were performed in R version 4.1.2.

Spatial and Temporal Variation in Salinity
* Salinity data were analyzed with ?

Multivariate Analyses

Community data from the transect surveys and field exclusion experiment were analyzed using the vegan package version 2.5-7. Species abundances were relativized with a double Wisconsin transformation. This standardizes species to equal maxima, then sites to equal totals, putting equalize emphasis among sample units and among species. We ordinated the data with non-metric multidimensional scaling (nMDS) plots. We then performed a permutational multivariate analysis of variance (PERMANOVA) to test the null hypothesis that the centroids of the groups are equivalent. In order to respect the dependence of sampling blocks within a site, the permutation scheme was restricted such that all quadrats along a transect, and plots within a site, for the surveys and field experiment respectively, were always permuted together. In doing so, 199 permutations were run on a matrix of Bray-Curtis dissimilarities. Because PERMANOVA cannot distinguish between differences in centroid location or levels of dispersion, we also used a PERMDISP test with 999 permutations, to test if the variances of the groups are different.
- Field exclusion experiment: Inclusion plots were not effective at containing L. pelta in the low salinity sites, therefore we only analyzed exclusion and control plots. A heat wave took place in August 2011 resulting in the die-off of many species, thus we only analyzed communities in the penultimate sampling point of July.
We plotted community trajectory data using the ecotraj package, version 0.0.1, which reflect the mean of the centroid for each of the 6 sites at each time point.

Univariate Analyses

Transect Surveys
Species abundance over time in low and high salinity sites were modelled with generalised linear models using the glmmTMB package version 1.1.2.3. Each species was modelled with month and salinity region as fixed factors, and a random effect of site nested within salinity region. For mobile invertebrates that were counted as total individuals (compared to percent cover), models were fit with with a negative binomial and poisson error distribution, with and without zero inflation. Sessile invertebrates and algal species were counted as percent cover of quadrats, and therefore we also fit models with a beta and a tweedie error distribution to these species. Models were then compared with the Akaike information criterion (AIC), and whichever had the lowest AIC value was chosen. Model diagnostics (residuals plots and Durbin-Watson tests) were run using the DHARMa package, version 0.4.4 (Hartig 2020). If variance across groups was unequal, we included a dispersion formula for either month, salinity region or both. Limpets, littorines, Chthamalus dalli, Red algae, Green algae and Brown algae were all fitted with a negative binomial distribution. Balanus glandula was fitted with a beta error distribution. To deal with the excessive number of zeros in the limpet count data, we also modelled these data with a zero-inflated negative binomial model. The significance of fixed effects was then tested using Wald $X^{2} tests in the car package version 3.0-11. Differences among the levels of the effects were tested using post-hoc pairwise comparisons, based on estimated marginal means using a Tukey adjustment with the emmeans package, version 1.6.3.

Field Exclusion Experiment
To test whether our exclusion plots were effective at keeping grazers out, we pooled all grazer species (limpets and littorines) and fit a generalized linear model with a negative binomial error distribution and a dispersion formula for salinity region. Because all months prior to July were likely to impact community structure, we included data from May, June and July in this model only. We fit generalized linear models to B. glandula, C. dalli, all algae pooled together, algae split into taxonomic groups of red, brown and green, as well as split into species which are grazed and non-grazed by limpets. There was minimal amounts of red and brown algae, thus we are only including models for green algae, and all algal species pooled together. C. dalli, all algae species, and green algae were all modelled as zero-inflated data as well. We were unable to fit a model including a random intercept of site nested within region to the B. glandula data without resulting in convergence issues. To resolve this, we fit the model with site as a fixed effect instead.

Results

Spatial and Temporal Variation in Salinity
Surface salinities showed a seasonal pattern across all sites, with highest salinities occurring in the winter and early spring months (October to April) and lowest salinities occurring in the summer and early autumn (June to September) (Fig. 2). Salinity in the Gulf Islands was consistently higher than in West Vancouver (ANOVA, P<0.0001), and this difference was greatest in the summer months (May to August).

Figure 2: Measured surface salinity (psu) from sites in the Gulf Islands (shaded) and in West Vancouver (unshaded), British Columbia. Dashed line indicates Fraser River discharge rate (10³m³/s) measured at Hope, British Columbia (Environment Canada, 2012). Surface salinity for Eagle Cove, April 7, 2011, was influenced by heavy rainfall.

Multivariate Analyses

Transect Surveys

A stable NMDS plot with 3 dimensions and a stress of 0.15 was reached after 30 tries. PERMANOVA showed a significant effect of salinity region (F = 81.176, p = 0.005; Table 1), month (F = 4.481, p = 0.005; Table 1), and the interaction between the two (F = 1.697, p = 0.005; Table 1) on community composition.

Plots from low and high salinity sites closely cluster together. In the low salinity sites, surveys from July and August begin to form slightly less variable clusters, which may be explained by the seasonal change in salinity at these sites.

Figure 3: nMDS plot of community data collected from transect surveys from 3 high and 3 low salinity sites, collected summer 2011

Table 1: Results of PERMANOVA analysis on Bray-Curtis dissimilarities for salinity region and month effects on intertidal communities.
treatments	df	SS	MS	pseudo.F	R2	p.value
salinity region	1	14.76	14.76	81.18	0.29	0.005
month	3	2.44	0.81	4.48	0.05	0.005
salinity region x month	3	0.93	0.31	1.7	0.02	0.005
residuals	184	33.46	0.18		0.65
total	191	51.60			1.00

A PERMDISP test revealed that dispersion among salinity groups was equal (F = 0.23, P > 0.05; table 2), but that May and June have less variance compared to September (F = 6.515, p = 0.001; table 3). As groups with different dispersions may result in misleadingly low p values, PERMANOVA results must be interpreted cautiously.

Figure 4: dispersion levels of centroids among salinity levels and during each month of sampling

Table 2: Results of PERMDISP analysis on Bray-Curtis dissimilarities for salinity region
Df	SS	MS	F	N.Perm	Pr..F.
1	0.0041347	0.0041347	0.229678703154221	999	0.662
190	3.4204103	0.0180022

Table 3: Results of PERMDISP analysis on Bray-Curtis dissimilarities for month
Df	SS	MS	F	N.Perm	Pr..F.
3	0.1408684	0.0469561	6.51498911371278	999	0.002
188	1.3549909	0.0072074

Simper analysis
Which species contribute at least 70% of the differences between groups?

Table 4: Results of Simper Analysis. Species that contribute the most to Bray-Curtis dissimilarities for salinity regions
species	contribution	p_value	group
mytilus_pt	0.1590332	0.115	salinity
balanus_pt	0.0862986	0.495	salinity
fucus_pt	0.0848790	0.115	salinity
chthamalus_pt	0.0811211	0.115	salinity
mastocarpus_crust_pt	0.0434617	0.585	salinity

Table 5: Results of Simper Analysis. Species that contribute the most to Bray-Curtis dissimilarities for different months
species	contribution	group	p_value
mytilus_pt	0.1329027	July_September	0.025
balanus_pt	0.1067202	July_September	0.005
fucus_pt	0.0835693	July_September	0.035
chthamalus_pt	0.0544253	July_September	0.915
ulva_pt	0.0464845	July_September	0.005
mytilus_pt	0.1062865	June_July	0.990
balanus_pt	0.0907696	June_July	0.145
fucus_pt	0.0827536	June_July	0.070
chthamalus_pt	0.0687799	June_July	0.105
barnacle_recruits_pt	0.0471294	June_July	0.105
mytilus_pt	0.1283783	June_September	0.085
balanus_pt	0.0867641	June_September	0.340
fucus_pt	0.0797205	June_September	0.220
chthamalus_pt	0.0698817	June_September	0.055
barnacle_recruits_pt	0.0471294	June_September	0.105
mytilus_pt	0.1127107	May_July	0.845
balanus_pt	0.0945008	May_July	0.035
fucus_pt	0.0774119	May_July	0.455
chthamalus_pt	0.0531743	May_July	0.905
pelta_no	0.0453223	May_July	0.015
mytilus_pt	0.1068071	May_June	0.970
barnacle_recruits_pt	0.0677211	May_June	0.005
chthamalus_pt	0.0671705	May_June	0.080
fucus_pt	0.0658586	May_June	0.995
balanus_pt	0.0550453	May_June	1.000
mytilus_pt	0.1335991	May_September	0.005
balanus_pt	0.0853340	May_September	0.530
fucus_pt	0.0741613	May_September	0.800
chthamalus_pt	0.0522939	May_September	0.895
ulva_pt	0.0454040	May_September	0.005

A critique of simper analysis is that generally the species with the most variance also have the highest contribution. Even if you make groups that are copies of each other, the method will single out species with high contribution, but these are not contributions to non-existing between-group differences but to random noise variation in species abundances. Here I’m plotting species variance across all plots against their calculated contributions (Figure 5).

Figure 5: variance vs contribution to dissimilarity matrix of each species

Community Trajectories

The three low salinity sites experience dips in salinity seasonally as the Fraser River discharge increases. If salinity is a major determinant in the structure of intertidal communities, either directly through physiological tolerance of species or indirectly through changes in species interactions, we would expect larger changes in the trajectories of the three sites located closer to the FRE. In order to visualize this, we ordinated each sites trajectories in nMDS space (Figure 6). The three high salinity sites are all clustered close together and only diverge slightly throughout summer, whereas the three low salinity sites vary considerably at each time point.

Figure 6: Trajectories of intertidal communities at low and high salinity sites during the summer of 2011. Species abundances were averaged across 8 quadrats along a transect at each site and time point, and a matrix of Bray-Curtis dissimilarities was consructed.

Field Exclusion Experiment
Species abundance: A stable NMDS plot with 3 dimensions and a stress of 0.11 was reached after 20 tries. PERMANOVA analysis revealed significant differences in the community composition of salinity region, grazer treatment and the interaction between the two.

We hypothesized that the exclusion plots in the high salinity sites would resemble the community structure in both the control and exclusion plots of the low salinity sites, due to the lack of limpet grazing. Both treatment (pseudo-F = , p value = 0.005; Table 6), as well as the interaction between salinity region and treatment (pseudo-F, p value = 0.025; Table 6) were significant.

Figure 7: nMDS plot of community data from field exclusion experiment in high and low salinity sites from the summer of 2011.

Table 6: Results of PERMANOVA analysis done on Bray-Curtis dissimilarities for salinity region and grazer manipulation effects on intertidal communities.
treatments	df	SS	MS	pseudo.F	R2	p.value
salinity region	1	0.97	0.97	8.2	0.09	0.375
treatment	1	0.40	0.4	3.41	0.04	0.005
salinity region x treatment	1	0.18	0.18	1.52	0.02	0.025
residuals	80	9.43	0.12		0.86
total	83	10.97			1.00

A PERMDISP test revealed that dispersion among salinity groups was unequal (F = 17.3, P = 0.001; table 6), as well as dispersion among treatment groups (F = 7.5, p = 0.007; table 7). As groups with different dispersions may result in misleadingly low p values, PERMANOVA results must be interpreted cautiously.

Figure 8: dispersion levels of centroids among salinity levels and treatment levels

Table 6: Results of PERMDISP analysis on Bray-Curtis dissimilarities for salinity region
Df	SS	MS	F	N.Perm	Pr..F.
1	0.4511003	0.4511003	17.3184440175554	999	0.001
82	2.1358863	0.0260474

Table 7: Results of PERMDISP analysis on Bray-Curtis dissimilarities for treatment
Df	SS	MS	F	N.Perm	Pr..F.
1	0.1553742	0.1553742	7.46326674014969	999	0.008
82	1.7071186	0.0208185

Simper analysis
How much does each species contribute to the overall dissimilarity matrix among salinity regions (Table 8), and treatments (Table 9).

Table 8: Results of Simper Analysis. Species that contribute the most to Bray-Curtis dissimilarities for salinity regions
species	contribution	p_value	group
ulva_pt	0.2125875	0.115	salinity
balanus_no	0.1133128	0.390	salinity
chthamalus_no	0.0857018	0.305	salinity
masto_crust_pt	0.0630440	0.685	salinity
fucus_pt	0.0556836	0.115	salinity

Table 9: Results of Simper Analysis. Species that contribute the most to Bray-Curtis dissimilarities for treatments
species	contribution	p_value	group
ulva_pt	0.1927748	0.005	treatment
balanus_no	0.1076078	0.010	treatment
chthamalus_no	0.0872806	0.005	treatment
masto_crust_pt	0.0622233	0.785	treatment
fucus_pt	0.0546080	0.095	treatment

Univariate Analyses

Transect Surveys

Invertebrates
All limpet species were pooled and include primarily Lottia paradigitalis, but also Lottia digitalis, Lottia pelta, Lottia scutum, and Lottia persona. Limpet abundance is consistently higher in the high salinity region, and saw a decreased over the course of the summer in the low salinity region only. Littorines were also pooled together and consisted of Littorina sitkana and Littorina scutulata. Littorine abundance was only higher in the high salinity region during the month of July, before decreasing slightly; Littorines also experienced a small dip in abundance during July in the low salinity region. Abundance of Balanus glandula increased over the summer in both regions, and Chthamalus dalli was consistently higher in the high salinity sites (Figure 9).

Figure 9: Invertebrate species surveyed in 3 high and 3 low salinity sites during summer 2011

Table 10: ANOVA test results, invertebrates from transects surveys
source	statistic	df	p.value
limpets
month	22.2133596	3	0.0000589
salinity_regime	13.2517681	1	0.0002723
month:salinity_regime	27.5102771	3	0.0000046
littorines
month	26.7784889	3	0.0000066
salinity_regime	0.0009192	1	0.9758133
month:salinity_regime	27.9464824	3	0.0000037
balanus
month	41.1194879	3	0.0000000
salinity_regime	1.2372062	1	0.2660102
chthamalus
month	4.9897982	3	0.1725457
salinity_regime	66.5440393	1	0.0000000

Algae by taxonomic group

Red Algae consists of primarily Pyropia, and Mastocarpus (foliose and crust form), but also Hildenbrandia, Endocladia, Microcladia, and Polysiphonia. Green algae consists of primarily Ulva spp, but also Urospora, and Acrosiphonia. Brown algae consists primarly of Fucus distichus, but also Melanosiphon spp, and Leathesia spp. There is weak evidence of a decrease in red algae over the summer in both regions (F= 6.26, p = 0.09; table 11), but was consistently higher in the high salinity region (F = 13.25, p < 0.001; table 11). Green and brown algae increased in both regions over the summer (Figure 10).

Figure 10: Algae divided into taxonomic groups, surveyed at 3 high and 3 low salinity sites during the summer of 2011.

Table 11: ANOVA test results, algae by taxonomy from transects surveys
source	statistic	df	p.value
brown algae
month	15.8875372	3	0.0011958
salinity_regime	1.3204216	1	0.2505164
green algae
month	37.8466010	3	0.0000000
salinity_regime	0.4019598	1	0.5260789
red algae
month	6.2589556	3	0.0996691
salinity_regime	13.2491129	1	0.0002727

Field Exclusion Experiment

Grazers
Grazer species are Lottia paradigitalis, Lottia pelta, Lottia digitalis, Lottia scutum, Littorina scutulata, and Littorina sitkana. Grazers were always more abundant in the high salinity sites, but decreased in numbers in the the exclusion plots (fig. 11).

figure 11: Grazer abundance in grazer exclusion vs. control plots in high and low salinity sites.

Table 11: ANOVA test results, all grazers from field exclusion experiment
source	statistic	df	p.value
treatment	0.0612935	1	0.8044631
region	10.0465690	1	0.0015263
treatment:region	7.7226311	1	0.0054533

Barnacles

There were site level differences in B. glandula (F = 65.31, p < 0.001; Fig.12), and weak evidence of a treatment level effect (F = 2.64, p = 0.10; Fig.12) . C. dalli decreased in the grazer exclusion plots in the high salinity sites only.

Figure 12: Barnacle abundance in grazer exclusion vs. control plots in high and low salinity sites.

Table 13: ANOVA test results, balanus glandula from field exclusion experiement
source	statistic	df	p.value
treatment	2.640882	1	0.1041457
site	65.318524	5	0.0000000
treatment:site	8.834319	5	0.1158579

Table 14: ANOVA test results, Chthamalus dalli from field exclusion experiement
source	statistic	df	p.value
treatment	0.0477698	1	0.8269905
region	2.7900916	1	0.0948488
treatment:region	17.1414786	1	0.0000347

All algae

Algae percent cover increased in the grazer exclusion plots in the high salinity sites (figure 13). This includes reds, greens, browns and diatoms.

Figure 13: Algae abundance in grazer exclusion vs. control plots in high and low salinity sites.

Table 15: ANOVA test results, All algae species pooled together, from field exclusion experiement
source	statistic	df	p.value
treatment	1.243033	1	0.2648873
region	19.010591	1	0.0000130
treatment:region	19.127575	1	0.0000122

Green Algae
Green algae in the experimental plots consists exclusively of Ulva. In the high salinity region, when grazers were excluded green algal cover increased by 32% (Fig. 14).

Figure 14: Green algae abundance in grazer exclusion vs. control plots in high and low salinity sites.

Table 16: ANOVA test results, All green algae species pooled together, from field exclusion experiement
term	statistic	df	p.value
treatment	0.6632273	1	0.4154229
region	18.6691845	1	0.0000155
treatment:region	17.5541920	1	0.0000279

Supplementary Material

Field Exclusion - multivariate

Because plots were cleared in May, species were unlikely to have reached carrying capacity by July, and species abundance may not accurately depict the impact of the experimental manipulation. To account for this in our analysis, we also collapsed the data to presence/absence and re-ran the ordination and permutation as described above, on a matrix of Jaccard dissimilarities.

Figure S1: nMDS plot of species richness data from field exlcusion experiment in high and low salinity sites.

A PERMANOVA on community composition of the presence or absence of species found treatment (pseudo F = 8.97, p value = 0.005; Table S1), salinity region (pseudo F = 19.8, p value = 0.005; Table S1) and the interaction between the two (pseudo F = 8.38, p value = 0.005; Table S1) were significant.

Table S1: Results of PERMANOVA analysis done on Jaccard dissimilarities presence and absence of species, for salinity region and grazer manipulation effects on intertidal communities.
treatments	df	SS	MS	pseudo.F	R2	p.value
salinity region	1	2.28	2.28	19.8	0.17	0.005
treatment	1	1.03	1.03	8.97	0.08	0.005
salinity region x treatment	1	0.97	0.97	8.38	0.07	0.005
residuals	80	9.22	0.12		0.68
total	83	13.51			1.00

Transect Surveys - univariate

Algae by functional group

We also grouped algae by functional groups. ‘Grazed’ algae are filaments and sheets that are more likely to be consumed by invertebrates at a higher frequency, and consist of Ulva, Pyropia, Polysiphonia, Diatoms, Urospora, Acrosiphonia. ‘Non-grazed’ algae are turfs and crusts that are consumed at a lower frequency, and consist of Fucus, Mastocarpus, Mastocarpus crust, Hildenbrandia, Endocladia, Melanosiphon, Leathesia. Grazed algae increased over the summer in the low salinity region only, and non-grazed algae increased over the summer in both regions (Figure S2).

Figure S2: Algae divided into functional groups, surveyed at 3 high and 3 low salinity sites during the summer of 2011.

Table S2: ANOVA test results, algae by functional group, from transect surveys
source	statistic	df	p.value
grazed algae
month	18.0978993	3	0.0004199
salinity_regime	2.5636368	1	0.1093468
month:salinity_regime	9.8419923	3	0.0199581
non grazed algae
month	17.2886876	3	0.0006164
salinity_regime	0.7150894	1	0.3977594
month:salinity_regime	3.4599121	3	0.3259987

Field Exclusion - univariate: